The crystallization of hard disks induced by a temperature gradient 
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While uniform temperature has no effect on equilibrium properties of hard-core systems, its 
gradient might substantially change their behaviour. In particular, in hard-disk system subject to 
temperature difference AT disks are repelled from the hot boundary of the system and accumulate 
at the cold one. Using event-driven molecular dynamics simulations we show that for sufficiently 
large AT or coverage ratio p* , crystal forms at the cold boundary. In this spatially inhomogeneous 
system a significant decrease of diffusivity of disks clearly marks the stationary interface between 
liquid and crystal. Such a behaviour is also supported through calculation of the radial distribution 
function and the bond order parameter. Simulations show that for this nonequilibrium system the 
equipartition of energy holds and velocity obeys the Boltzmann distribution. 



I. INTRODUCTION 



Hard-core systems, whose studies were initiated more 
than fifty years ago with a seminal papers of Alder and 
Wainwright [l|, l2|, are still interesting. The apparent 
simplicity of these models is very deceptive and their be- 
haviour even with respect to some equilibrium properties 
is not yet fully understood. An example of that is the 
nature of the liquid-crystal phase transition in the two- 
dimensional hard-disk system. Although the initial pro- 
posal that this transition should be discontinuous 3] got 
some support from density- functional approaches |^-u|, 
alternative explanation challenged this possibility. In the 
so-called KTHNY scenario [8l-tl3| the crystal melts into a 
liquid via an intervening " hexatic" phase and as a result 
the transition consists of two continuous transitions. Re- 
cent extensive simulations support the KTHNY scenario 
[lil UH and show Van der Waals loops, that would indi- 
cate a discontinuous transition are only finite size effects. 
Although various experimental works also agree with the 
KTHNY scenario |16h21| , the problem is not yet entirely 
settled 

Nonequilibrium behaviour is even more challenging es- 
pecially in low dimensional systems where transport coef- 
ficients are known to diverge and the validity of the con- 
ventional hydrodynamic description is at risk [23]. Such 
divergences are related with the rate of decay of some 
correlation functions. The first estimations by Alder in 
Wainwright in 1970, suggested their slow, power-law de- 
cay [24|. but it was initially considered with some scep- 
ticism 25]. However, subsequent simulations seemed to 
support the slow decay and recent extensive calculations 
modify Alder and Wainwright estimation only by some 
logarithmic corrections [26]. 

Nonequilibrium conditions are often induced through 
a temperature gradient and low-dimensional hard core 
systems were also analysed under such conditions. In 
particular, the behaviour of heat conductivity in one- 
dimensional systems has attracted considerable attention 
[27M29| | since novel, and to some extent universal, power 
laws were found to describe its divergence with the sys- 



tem size. As another example, one can mention spatial 
segregation, known as the Ludwig-Soret effect, that takes 
place in non mono-dispersed systems subject to temper- 
ature gradient in [3fJ, |31( . 

Hard-core systems provide fundamental paradigms for 
understanding liquids, gases, glasses, crystals, colloids 
and even colloidal crystals. Especially the latter ones are 
recently of much interest. This is because the fine control 
of their properties and behaviour is now possible due to 
novel techniques that use for example microgravity, di- 
electrophoresis, or colloidal epitaxy [32[ . Moreover, col- 
loidal crystals obtained with such techniques might find 
numerous applications as photonic crystals (33j , optic fil- 
ters [lij or chemical sensors [35j . 

Yet another technique used for fabrication of colloidal 
crystals relies on the use of the temperature gradient 
36]. This intrinsically nonequilibrium technique offers 
fine control of the growth process and fabricated in such 
away crystals are relatively large (size ~3mm) [37j ■ It 
is perhaps surprising that the temperature gradient in 
hard-core systems can induce crystal growth since uni- 
form temperature is in a sense irrelevant and does not 
determine their phase diagram (energy of interactions 
in hard-core systems is either zero or infinity and their 
transformations are entirely of an entropic origin). Such 
a counterintuitive but rather simple idea underlying the 
temperature-gradient method encourages numerical cal- 
culations that would support experimental findings. Such 
a corroboration might lead to a better understanding of 
the crystal growth that even in hard-core systems still 
offers some surprises [38[ . In the present paper we report 
results of event-driven simulations of crystallization in a 
system of hard-disks subject to a temperature gradient. 
In Sect. II we introduce the model and method. Results 
are described in Sect. Ill and we conclude in Sect. IV. 



II. MODEL AND METHOD 



In our model N identical hard disks of radius r and 
unit mass are moving in a box of size L x and L y . At 
the left and right ends of the channel there are thermal 
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walls that are kept at temperatures T\ and T9, respec- 
tively (Ti > T 2 ). In the thermostat we used [39|], after 
the collision with the wall kept at temperature T a disk 
has its normal component sampled from the distribution 

p( v x) = ^^g^s exp(— with the sign in the argument 
of the Heaviside function depending on the location of the 
wall. Its parallel component is sampled from the Gaus- 

v 2 

sian distribution p(v y ) = ^^j ex P(~2r)" ^ n ^ ne ver tical 
(y) direction periodic boundary conditions are imposed. 

An efficient way to simulate hard-disk systems is to use 
event-driven molecular dynamics [40l - l42T | . Performance of 
the algorithm considerably increases upon implementing 
heap searching and sectorization and such methods have 
already been applied to a number of problems [H, SUSHI- 

Initially, disks are uniformly distributed inside the 
box. Their velocity obeys the Maxwell-Boltzmann dis- 
tribution with spatially dependent temperature (i.e., ki- 
netic energy) interpolating linearly between T\ and T2 
(such a choice speeds up the relaxation of the systems). 
Usually square geometry was used (L x = L y ) and the 
number of disks N ranges from 900 (30x30) to 3600 
(60x60). In some cases rectangular geometry was used 
with N = 2700 (90x30) disks. Reported results do not 
depend significantly on the finite size effects and such 
modest number of disks, in our opinion, is sufficient to 
analyse crystallization induced by the temperature gra- 
dient, as described in the present paper. 

The behaviour of hard-disks systems is controlled by 
the coverage ratio defined as 



where Nirr 2 is the total area covered by disks and L x L y 
is the area of the two-dimensional simulation box. In our 
simulations p* varied between 0.125 and 0.754. As we 
will show below, the temperature difference AT = T\—Ti 
plays also an important role in determining the behaviour 
in our model. 



III. RESULTS 

Qualitatively, the evolution of our model is illustrated 
in FigJI] Starting from the uniform distribution, disks are 
gradually expelled from the hot (left) boundary and move 
to the cold (right) boundary. Sufficient accumulation of 
disks induces crystallization (hexagonal structure). Only 
a fraction of disks crystallize and the hot end remains 
disordered. Sometimes the resulting crystal is not perfect 
and some vacancies form that might persist for a very 
long time. 

Further evidence of crystal structure formation is seen 
in Fig|3]where the spatial dependence of the average local 
density of disks is presented. Pronounced oscillations at 
the right boundary indicate formation of crystal layers. 

To examine the behaviour of our model more quanti- 
tatively we measured the mean squared displacement of 
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FIG. 1. Time evolution of 2700 hard disks of density p* — 

0. 607 simulated in 90 x 30 box and for Ti = 5 and T 2 = 

1. Approximately at t = 10 at the right (cold) boundary 
crystallization starts. Around t — 50 the system reaches the 
stationary configuration. The single-site vacancy that is close 
to the right boundary at t = 10 is still visible also for t = 50. 

the disks (r 2 ). Since the system is spatially inhomoge- 
neous, we divided the simulation box into several slices 
and calculated (r 2 ) and other quantities for disks in each 
slice separately (disks contribute to a given slice if they 
start in this slice). To enable comparison of diffusive 
properties, we fixed the simulation time as t — 10. Such 
a value is much larger than mean time of ballistic (free) 
motion and usually smaller than the time where diffusion 
of disks becomes affected by finite size effects. 

In the following we discuss dependence of the diffu- 
sive behaviour on the coverage ratio p* , the tempera- 
ture difference AT, and size effects. To further analyse 
the local ordering in our system we measured the ra- 
dial distribution function and the bond order parameter. 
Finally, we show that this strongly nonequilibrium and 
inhomogeneous system locally can be treated as a quasi- 
cquilibrium one, where equipartition of energy holds and 
velocity of the disks obeys the Boltzmann distribution. 



A. Diffusion 

In Fig. [3] we present the mean squared displacement 
as a function of the distance L from the hot boundary. 
The system was divided into 30 slices and we measured 
separately x (horizontal) and y (vertical) components of 



3 




30 
20 
10 




20 40 60 80 
L 



t=10 



30 
20 
10 




20 40 60 80 
L 



t=50 



20 40 60 80 
L 



20 40 60 80 
L 



FIG. 2. Time evolution of the spatially dependent local den- 
sity. Simulations were made for the same system as in FigfJJ 
and the results are averaged over 100 independent runs. 



the displacement. Starting from the initial configuration 
the system was first relaxed for a time sufficient to reach 
the stationary state and after that measurement of the 
displacement was performed. 

These data show that for sufficiently large density, dif- 
fusion of disks close to the cold boundary drops to very 
small values. These are the disks that form the crystal. 
The remaining disks diffuse freely but their diffusivity 
depends on their position. The closer to the hot bound- 
ary the disks are, the larger their diffusivity is. However, 
very close to the hot boundary this tendency reveres. 
Such a reduction of diffusivity is a natural consequence 
of a proximity to the boundary. It is also interesting to 
notice that although our model is anisotropic, diffusion 
in the x and y directions is similar. 

The plot of displacement (r 2 ) on the logarithmic scale 
(FigEl bottom panel) enabled us to approximately lo- 
cate the onset of crystal formation. Namely, p* = 0.607 
seems to be a threshold value that separate data that at 
the cold boundary show negligible (~ 0.1) displacement 
from those with a larger value of (r 2 ). We will use such 
a criterion for the construction of the phase diagram of 
the model (FigJSJ). Let us also notice that the threshold 
value p* — 0.607 is much lower than the recently reported 
values (0.723) for hard disks without a temperature dif- 
ference (l4T |. 

That for hard disks and sufficiently large coverage ra- 
tio p* crystallization takes place is not surprising. Less 
obvious seems to be crystalization induced by a temper- 
ature difference. In Fig. 2] we show the displacement 
(r 2 ) calculated for the fixed coverage ratio (p* = 0.636) 
but with varying temperature difference AT. At the cold 
boundary a transition around T\ = 3 can be seen between 
low-Ti (concave) and high-Ti (convex) regimes that we 
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FIG. 3. Mean-squared displacement after time t = 10 as a 
function of distance from hot boundary for TV = 900, L x = 
L y — 30 and T\ = 5 and T2 = 1: (a) x component (a; 2 ) 
calculated for (curves from top to bottom): p* =0.581, 0.607, 
0.636, 0.664, 0.694, 0.723 and 0.754; (b) the y component 
(y 2 ) for the same p* as in (a); (c) the isotropic displacement 
(r 2 ) for the same p* as in (a); (d) the same as (c) but (r 2 ) is 
plotted on the logarithmic scale. Top four curves are obtained 
for (from top) p* = 0.125, 0.282, 0.502, and 0.554. Remaining 
curves are obtained for the same values as in (a) but with 
omitted value p* — 0.636. 



again identify as the formation of a crystal. 

Let us also notice that data in Fig|4] for different T\ 
intersect almost at the same point (L = 19). We do 
not understand such a behaviour yet, but assuming that 
for increasing T\ the displacement (r 2 ) will tend to the 
step-like function with a discontinuity at L = 19, we can 
interpret this point as a high-Ti limit of the existence 
of the crystal phase. What is in our opinion interesting, 
is the fact that location of this point can be estimated 
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FIG. 4. The mean-square displacement (r 2 ) after time t = 10 
as a function of distance from the hot boundary L calculated 
for p* = 0.636, N = 900, T 2 = 1, and T x equal to 1.1 (upward 
triangles), 1.5 (downward triangles), 2 (diamonds), 3 (circles), 
4 (squares), 5 (pentagons) and 10 (crosses). 



already from low-Ti data. It is also possible that the 
existence of such a point is related to the already men- 
tioned KTHNY scenario. Namely, moving from the cold 
to the hot boundary we effectively move from a crystal 
to a liquid phase. Our data e.g., for T\ = 10 show that 
diffusion vanishes for L — 24 and that marks the limit 
of existence of the crystal phase. Provided that the liq- 
uid phase terminates at L = 19 we can speculate that 
for 19 < L < 24 the system remains in the hexatic-likc 
phase. Certainly, to establish if this is true or not further 
analysis would be needed. 

Since the number of disks N in simulated systems is 
rather modest it is important to examine the role of 
finite-size effects. In Figl5]we present the displacement 
(r 2 ) calculated for several values of N. For fixed p* the 
change of N implies a change of the size L x (— L y ) and to 
enable comparison the horizontal axis in Fig(5]is rescaled. 
Only in the vicinity of the hot boundary there are some 
changes with the system size and in the rest of the system 
the simulation data already with N = 900 disks seems 
to be sufficient. Let us notice that for fixed tempera- 
ture difference (Ti — T2) and increasing system size, the 
temperature gradient decreases. One might expect that 
for decreasing temperature gradient the fraction of the 
system that crystallizes will also decrease. Simulations 
show that, at least for the examined systems, this is not 
the case. Our results (Figj5j) show fast convergence (with 
the system size N) and the resulting crystal occupies a 
finite fraction of the system. On the other hand, simu- 
lations with fixed temperature gradient (not presented) 
show that for increasing system size the size of the crystal 
also increases. Further analysis of the size dependence in 
hard disks systems with a temperature gradient would 
be desirable. 
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FIG. 5. The mean-square displacement (r 2 ) after time t = 10 
as a function of the reduced distance from the hot bound- 
ary L/L x calculated for p* — 0.636 (a) and p* = 0.694 (b). 
Simulations were made for Ti = 5, Ti = 1 and the number 
of disks N — 900 (upward triangles), N = 1600 (downward 
triangles), N = 2500 (diamonds), and N = 3600 (circles). 
Square geometry was used (L x — L y ). Let us notice that 
these data suggest the existence of a well-defined and station- 
ary^) boundary between crystal and liquid phases. 



B. Local ordering 

Substantial change of the diffusive behaviour could in- 
dicate crystal formation. In principle however, a drop 
of diffusivity might as well be due to formation of e.g., 
glassy phase. In the present subsection we analyse the 
local ordering in more details and that provide further 
evidence of the crystal formation. 

First, we show the results of the calculation of the ra- 
dial distribution function (Fig|6]). At the hot part of 
the system [L = Q.1L X ) only small oscillations are seen. 
They are more pronounced at the center (L = 0.51^) and 
become very strong at the cold part (L — 0.9La;). In the 
latter case the split of the second peak is seen and there 
are some hints that its formation might be considered as 
a precursor of the long-range order formation (fioT]). 

Further evidence of the crystal formation comes from 
the calculation of the bond order parameter [3143], (qe), 
where the average is taken over all disks in a given slab 
(for a more detailed definition of (qe) see e.g., (48|). In 
Figj7]we present (q 6 ) as a function of the reduced distance 
from the hot boundary L/L x calculated for p* — 0.636 
(upward and downward triangles) and p* — 0.694 (dia- 
monds and circles). Simulations were made for Ti = 5, 
T-2 = 1 and the number of disks N = 900 (upward trian- 
gles and diamonds) and TV = 1600 (downward triangles 
and circles). From these data one can see that in the hot 
part (qe) is close to zero and that is the expected value 
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FIG. 6. The radial distribution function as a function of the 
scaled distance R/r calculated for three slabs centered around 
L — 0.1L X (hot part), L — 0.5L X (center), and L = 0.9L X 
(cold part). Calculations were made for p* = 0.636, Ti = 5, 
T 2 = 1, and N = 3600 disks. 
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FIG. 7. Mean bond order parameter, (qe), as a function of the 
reduced distance from the hot boundary L/L x calculated for 
p* — 0.636 (upward and downward triangles) and p* = 0.694 
(diamonds and circles). Simulations were made for Ti = 5, 
Ti = 1 and the number of disks N = 900 (upward triangles 
and diamonds) and iV = 1600 (downward triangles and cir- 
cles). Square geometry was used (L x = L y ). The average is 
taken over all disks in a given sector. 



in the liquid phase. In the cold part (qe) is much larger 
which is yet another indicator of the formation of crystal 
structure. 

Let us finally point out that the region with positive 
diffusion (liquid) seems to be spatially well separated 
from the region with vanishing diffusion (crystal). Thus 
our model offers the possibility to study stationary prop- 
erties of crystal-liquid interface, however, more detailed 
analysis of this issue is yet to be done. 



C. Phase diagram 

Our main finding is that in a hard-disk system suffi- 
ciently strong temperature difference AT induces crystal 
formation. The threshold value of AT depends on p* and 
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FIG. 8. Phase diagram for the hard disk system of cover- 
age ratio p* and with temperature difference AT (T = 1). 
Simulations made for the number of disks N equal to 900 
(squares), 1600 (traingles), and 3600 (circles) gave nearly the 
same results. The onset of crystallization is detected on the 
value of displacement (r 2 ) for disks in the 10% of the system 
that is closest to the cold boundary. For a vanishing temper- 
ature difference (AT = 0) our threshold value p* « 0.72, is 
in a good agreement with the crystallization threshold 0.723 
estimated by Kozak et al. [22]]. 



the corresponding phase diagram is shown in FigjS] As 
a criterion for the appearance of the crystal we used the 
behaviour of the displacement (r 2 ). We considered that 
crystal was formed if disks in the 10% of the system that 
is closest to the cold boundary had a displacement (r 2 ) 
(after time t = 10) below 0.1. A noticeable feature is the 
negligible role of finite size effects and accurate location 
of the liquid-crystal boundary can be made already with 
simulations of N = 900 disks. 

One of the open problems related with this phase dia- 
gram concerns a possible existence of the hexatic phase 
in the AT > regime, and we already mentioned such a 
possibility discussing Fig|4] 



D. Quasi-equilibrium 

Although nonequilibrium systems posses a number of 
distinctive features, recent studies show that some no- 
tions and concepts that we use to describe equilibrium 
systems might be applicable also to nonequilibrium ones. 
Indeed, there is a number of works where nonequilibrium 
generalizations of entropy, temperature, or free-energy 
were successfully used. However, equilibrium analogies 
might be misleading as well. In particular, there are some 
indications that various, and in equilibrium equivalent, 
definitions of temperature in nonequilibrium yield dif- 
ferent values [4^ ]. Intensively studied systems with such 
ambiguities include driven granular mixtures, where each 
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FIG. 9. (a) Average square of horizontal (downward trian- 
gles) and vertical (upward triangles) components of velocity 
as a function of the distance from the hot boundary L. Sim- 
ulations are made for N = 900 disks, Ti = 5 and T% = 1 and 
for coverage ratio p* = 0.636. The system is divided into 30 
sectors and each point is an average for disks in a given sector. 
Within numerical accuracy these data support equipartition 
of energy in this model; (b) Distribution of \v x | in four choosen 
sectors located at L = 2 (upward triangles), L = 10) (down- 
ward triangles), L = 20) (diamonds), and L = 28) (circles). 
The inset shows the distribution in the semi-logarithmic plot 
and the linearity confirms the Boltzmann distribution. 



component seems to be characterized by its own tempera- 
ture [50( | . Moreover, in some glassy systems the effective 
temperature defined through fluctuation-dissipation re- 
lations turns out to depend on the choice of observables 
_51j. Some other cornerstones of equilibrium statistical 
mechanics like equipartition of energy (52j or gaussian 
distribution of velocities f53j were also reported to break 
down in some granular systems. 

The nonequilibrium character of our model is caused 
by the temperature gradient imposed at the boundaries. 
Nevertheless, it is tempting to think of such a system as 
being in local equilibrium. In this subsection we show 
that this is indeed the case, at least with respect to the 
equipartition of energy and distribution of velocities. In 
the upper part of FiglH] we present the average square 
of the horizontal ((v^)) and vertical ((«„)) components 
of the velocity. According to equipartition of energy in 



equilibrium systems they should be equal 

2 .2 2 ' . 

Since our model is nonequilibrium and spatially 
anisotropic (heat transfer takes place only in the hori- 
zontal direction) the equality (v x ) = (vy) is by no means 
obvious (let us recall that in the perpendicular (y) di- 
rection we imposed periodic boundary conditions). Our 
simulations, however, clearly support equipartition of en- 
ergy in this case (Fig [9]). Let us also notice that although 
for such a temperature difference and coverage ratio the 
crystal forms at the cold (right) boundary of the sys- 
tem, no singularity in the behaviour of (v x ), (v^) is seen 
(from Figf5] it follows that in this case crystal exists for 
(24 < L < 30). 

In the lower part of Fig[9] we present distribution of 
\v x \ in four sectors (only the distribution for v x is pre- 
sented and the results for v y are qualitatively similar). 
The familiar Boltzmann distribution is confirmed in the 
inset of Fig|9] where the distribution is plotted in the 
semi-logarithmic scale. When the sectors move from left 
to right (from hot to cold) both the variance of the dis- 
tribution and average \v x \ decrease. 



IV. CONCLUSIONS 

In the present paper, using event-driven molecular dy- 
namics simulations, we examined the hard disk system 
with coverage ratio p* and subject to temperature dif- 
ference AT. Our results show that in such a system the 
disks accumulate at the cold boundary and for sufficiently 
large p* or AT they crystallize at that side. We measured 
diffusion of disks along the system, density profiles, ra- 
dial distribution function, and bond order parameter and 
found that there is a clear separation of liquid and crys- 
tal. In hard-disk models, controlled only by coverage 
ratio p* , the system is either in the crystal or the liquid 
phase and their interface usually has only a transient na- 
ture. In our model such an interface is stationary and 
its further analysis might provide interesting results on 
the coexistence of liquid and crystal. However, further 
analysis of this issue we left for the future. Another in- 
teresting question left for the future concerns a possible 
existence of the hexatic phase in hard disks with the tem- 
perature gradient. If it exists, then again, in our model 
we would have a possibility to study stationary coexis- 
tence of this phase with the crystal and liquid phases. 
It would be also desirable to extend our work to three- 
dimensional hard-sphere system that would be closest to 
experimental realizations. 
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